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[i]  Reliable  predictions  of  groundwater  flow  and  solute  transport  require  an  estimation  of 
the  detailed  distribution  of  the  parameters  (e.g.,  hydraulic  conductivity,  effective  porosity) 
controlling  these  processes.  However,  such  parameters  are  difficult  to  estimate  because  of 
the  inaccessibility  and  complexity  of  the  subsurface.  In  this  regard,  developments  in 
parameter  estimation  techniques  and  investigations  of  field  experiments  are  still  challenging 
and  necessary  to  improve  our  understanding  and  the  prediction  of  hydrological  processes. 

Here  we  analyze  a  conservative  tracer  test  conducted  at  the  Boise  Hydrogeophysical 
Research  Site  in  2001  in  a  heterogeneous  unconfined  fluvial  aquifer.  Some  relevant 
characteristics  of  this  test  include:  variable-density  (sinking)  effects  because  of  the  injection 
concentration  of  the  bromide  tracer,  the  relatively  small  size  of  the  experiment,  and  the 
availability  of  various  sources  of  geophysical  and  hydrological  information.  The 
information  contained  in  this  experiment  is  evaluated  through  several  parameter  estimation 
approaches,  including  a  grid-search-based  strategy,  stochastic  simulation  of  hydrological 
property  distributions,  and  deterministic  inversion  using  regularization  and  pilot-point 
techniques.  Doing  this  allows  us  to  investigate  hydraulic  conductivity  and  effective  porosity 
distributions  and  to  compare  the  effects  of  assumptions  from  several  methods  and 
parameterizations.  Our  results  provide  new  insights  into  the  understanding  of  variable- 
density  transport  processes  and  the  hydrological  relevance  of  incorporating  various  sources 
of  information  in  parameter  estimation  approaches.  Among  others,  the  variable-density 
effect  and  the  effective  porosity  distribution,  as  well  as  their  coupling  with  the  hydraulic 
conductivity  structure,  are  seen  to  be  significant  in  the  transport  process.  The  results  also 
show  that  assumed  prior  information  can  strongly  influence  the  estimated  distributions  of 
hydrological  properties. 
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1.  Introduction 

[2]  Accurate  prediction  of  flow  and  solute  transport  in 
aquifers  is  a  key  prerequisite  for  effective  sustainable  man¬ 
agement  and  remediation  of  groundwater  resources.  In  fact, 
such  predictions  are  very  challenging  because  of  the  com¬ 
plexity  of  the  hydrological  processes  and  the  difficulty  in 
obtaining  information  about  important  flow  parameters  in 
the  subsurface.  The  general  procedure  for  investigating  the 
distribution  of  hydrological  properties,  such  as  hydraulic 
conductivity  ( K)  and  effective  porosity  (<5),  for  example, 
includes  (1)  experimental  design,  (2)  data  acquisition,  (3) 
model  conceptualization,  and  (4)  parameter  estimation 
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[e.g.,  de  Marsily  et  al.,  1999;  Carrera  et  al.,  2005].  First, 
designing  an  experiment  involves  the  development  of  data 
acquisition  techniques  and  strategies  that  are  sensitive  to 
properties  and  processes  of  interest.  Such  strategies  also 
need  to  be  optimized  with  regard  to  actual  field  and  bound¬ 
ary  conditions  as  well  as  cost  and  time  factors.  Second,  the 
objective  of  data  acquisition  is  to  gain  as  much  reliable  and 
useful  information  as  possible  for  the  goals  of  a  given 
study.  Third,  model  conceptualization  includes  identifying 
the  nature  of  the  problem  and  hydrogeologic  system,  select¬ 
ing  the  governing  equations,  defining  the  boundary  condi¬ 
tions  and  the  time  regime,  selecting  the  information  sources 
of  interest,  and  the  way  of  integrating  them.  This  step  is 
highly  dependent  on  the  prior  understanding  of  the  system 
and  processes.  Fourth  and  finally,  the  parameter  estimation 
is  focused  on  assigning  values  to  the  model  properties  iden¬ 
tified  as  controlling  the  processes  of  interest.  This  step  is 
often  performed  through  inverse  modeling  or  simulation 
approaches  to  infer  distributions  of  properties  that  allow 
matching  (to  an  acceptable  tolerance)  of  measurements  of 
the  processes  of  interest. 
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[3]  Advancements  in  the  above  areas  are  supported  by 
numerical,  laboratory,  and  field  experiments  [e.g.,  Rubin, 
1995;  Yeh  and  Liu,  2000 ;  Cardiff  et  al. ,  2009 ;  Pollock  and 
Cirpka,  2010;  Sudicky  et  al.,  2010],  In  particular,  field 
experiments  in  various  hydrological  and  geological  con¬ 
texts,  at  various  scales,  and  with  various  survey  strategies, 
are  essential  to  better  understand  flow  and  transport  proc¬ 
esses,  validate  and  improve  modeling  approaches,  and  guide 
novel  developments  [e.g.,  Sudicky,  1986;  Ptak  and  Schmid, 
1996;  Illman  et  al.,  2009;  Cvetkovic  et  al.,  2010;  Brauchler 
et  al.,  2011],  Such  field  experiments  commonly  include  var¬ 
iations  on  tracer  tests,  multiwell  pumping  or  injection  tests, 
and  single-well  slug  or  flowmeter  tests  or  direct-push  pro¬ 
files  [e.g.,  Hyndman  and  Gorelick,  1996;  Brauchler  et  al., 
2003;  Doherty,  2003;  Zhu  and  Yeh,  2005;  Cardiff  et  al., 
2011].  In  this  study  we  present  the  analysis  of  a  tracer  test 
performed  in  2001  at  the  Boise  Hydrogeophysical  Research 
Site  (BHRS)  [Barrash  et  al.,  2002;  Hausrath  et  al.,  2002], 
and  investigate  a  number  of  issues  related  to  the  transport  of 
solute  subject  to  variable-density  effects  in  a  heterogeneous 
unconfined  aquifer  at  a  local  scale  with  several  parameter 
estimation  approaches. 

[4]  One  main  objective  of  studies  about  solute  transport 
[e.g.,  Mackay  et  al.,  1986;  Sudicky,  1986;  Leblanc  et  al., 
1991 ;  Hyndman  et  al.,  2000;  Vereecken  et  al.,  2000;  Hub¬ 
bard  et  al.,  2001;  Scheibe  and  Chien,  2003;  Miiller  et  al., 
2010;  Sudicky  et  al.,  2010]  is  to  improve  the  estimation  of 
spatially  variable  properties  controlling  flow  and  transport, 
primarily  hydraulic  conductivity  (K)  and  effective  porosity 
(4>),  and  to  improve  the  understanding  and  prediction  of 
processes,  such  as  transport,  dispersion,  and  chemical  and 
biological  reactions  among  others  [e.g.,  de  Marsily  et  al., 
2005],  In  this  study,  we  evaluate  the  spatial  distribution  of 
hydrological  properties  (e.g.,  K,  </>,  dispersivity)  controlling 
the  observed  transport  process,  the  importance  of  each  one 
in  predicting  this  process,  and  the  possible  coupling  effects 
between  properties.  The  resulting  understanding  is  useful 
not  only  to  improve  predictions  or  solve  issues  at  one  spe¬ 
cific  site,  but  also  to  allow  the  extrapolation  of  gained 
knowledge  to  other  areas  where  only  limited  information 
may  be  available  [e.g.,  Hyndman  et  al.,  2000;  Linde  et  al., 
2006 ;  Dafflon  et  al.,  2010]. 

[5]  Another  subject  of  interest  in  this  study  is  how  we 
can  use  indirect  or  supplemental  information  to  constrain 
the  distribution  of  hydrological  properties  controlling  flow 
and  transport  processes  and  thus  improve  flow  and  trans¬ 
port  predictions.  In  this  context,  geophysical  methods  have 
shown  much  potential  [e.g.,  Hyndman  and  Gorelick,  1996; 
Slater  et  al.,  1997 ;  Hubbard  et  al.,  2001 ;  Singha  and  Gorelick, 
2006;  Muller  et  al.,  2010]  because  they  provide  a  scale  of 
spatial  resolution  and  degree  of  subsurface  coverage  not 
available  with  traditional  hydrological  measurement  techni¬ 
ques  (e.g.,  borehole  logs,  core  analyses,  pumping,  and  tracer 
tests).  In  addition,  geophysical  methods  are  able  to  provide 
large  quantities  of  informative  data  in  less  time  and  at  lower 
cost  than  most  hydrological  methods.  The  BHRS  is  particu¬ 
larly  suitable  for  investigating  the  value  of  supplemental  in¬ 
formation  due  to  the  broad  database  of  geophysical  and 
hydrological  information  available  at  this  site.  In  this  study 
we  examine  how  hydrological  information  gained  from  a 
tracer  test  alone  compares  with  information  inferred  from 
geophysical  data,  and  how  the  4>  distribution  estimated 


from  geophysical  data  can  influence  the  hydrological  pa¬ 
rameter  estimation  process.  This  is  essential  to  evaluate 
how  geophysical  data  may  be  sensitive  to  hydrological 
properties,  how  such  properties  may  be  related  together,  and 
what  can  be  expected  from  various  hydrological,  geophysi¬ 
cal,  and  hydrogeophysical  data  sets. 

[6]  In  addition,  the  use  of  parameter  estimation  techni¬ 
ques  for  variable-density  solute  transport  studies  are  still 
computationally  intensive  and  have  not  been  investigated 
much  with  regard  to  the  information  content  of  hydrologi¬ 
cal  data  sets  influenced  by  density-driven  processes. 
Although  a  number  of  studies  have  recognized  or  described 
sinking  of  a  tracer  because  of  density  effects  [e.g.,  Leblanc 
et  al.,  1991 ;  Vereecken  et  al.,  2000;  Beinhorn  et  al.,  2005 ; 
Miiller  et  al,  2010],  the  use  of  parameter  estimation  meth¬ 
ods  for  problems  that  include  density  effects  is  still  rare. 
Indeed,  it  has  been  noted  that  a  significant  number  of  stud¬ 
ies  neglect  the  variable-density  effect  although  they  should 
not  [ Simmons ,  2005].  This  is  probably  because  of  the 
increased  complexity  in  modeling  such  processes  and  in 
processing  and  inverting  such  data,  to  significantly  increased 
demand  for  computing  power,  and  to  increased  difficulty  in 
interpretation  of  the  results.  Here  for  inverting  such  data 
sets,  we  investigate  several  parameter  estimation  approaches, 
and  some  crucial  issues  such  as  the  optimization  of  comput¬ 
ing  time  and  model  size,  the  objective  function  used,  and  the 
implementation  of  inversion  parametric  constraints  (i.e., 
“prior”  information). 

[7]  This  paper  is  organized  as  follows.  First,  the  BHRS 
and  hydrogeological  and  geophysical  information  relevant 
for  the  tracer  test  are  presented.  Then  the  tracer  test  is 
described  as  well  as  the  parameters  used  for  the  flow  and 
solute  transport  modeling.  Next,  the  measured  concentra¬ 
tions  during  the  tracer  test  are  used  to  evaluate  the  distribu¬ 
tion  of  hydrological  properties  such  as  K  and  <j>,  and  also  to 
investigate  the  information  contained  in  such  measure¬ 
ments  through  various  strategies  including  (1)  a  grid- 
search-based  approach  to  infer  relatively  simple  models  of 
hydrological  properties,  (2)  simulations  of  hydrological 
property  distributions  to  consider  the  impact  of  introducing 
small-scale  heterogeneity  for  predicting  the  measured  con¬ 
centrations,  and  (3)  a  deterministic  inversion  involving  reg¬ 
ularization  and  pilot-point  techniques  to  evaluate  various 
parameterizations  of  the  problem. 

2.  Field  Site  Hydrogeology 

[8]  The  BHRS  is  a  hydrological  and  geophysical  field 
research  site  located  on  a  gravel  bar  adjacent  to  the  Boise 
River  ~15  km  from  downtown  Boise,  Idaho.  The  shallow, 
unconfined  aquifer  at  the  site  consists  of  coarse  unconsoli¬ 
dated  fluvial  deposits  that  are  ~20-m-thick,  have  minimal 
fractions  of  silt  and  clay,  and  are  underlain  by  a  layer  of 
red  clay  [ Barrash  and  Clemo,  2002 ;  Barrash  and  Reboulet, 
2004].  At  this  site,  18  wells  were  emplaced  carefully  using 
a  core  drive-drill  method  to  minimize  the  disturbance  of 
the  surrounding  formation  and  were  fully  screened  through 
the  aquifer  with  4-in  PVC.  The  well  field  consists  of  13 
wells  in  the  central  area  (~20  m  in  diameter)  and  five 
boundary  wells  ~  10-35  m  from  the  central  area  (Figure  1). 
The  arrangement  of  the  13  inner  wells  is  a  central  well 
(Al)  surrounded  by  two  concentric  rings  of  six  wells  each 
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Figure  1.  Aerial  view  of  the  BHRS  and  detailed  map  of 
the  central  area  wells  with  lines  to  indicate  between  which 
wells  the  tracer  test  (black,  solid)  and  static  cross-hole  GPR 
(gray)  acquisitions  have  been  performed.  The  time-lapse 
cross-hole  GPR  profile  shown  in  this  study  is  indicated 
with  a  dotted  line. 


(B1-B6  and  C1-C6).  The  distances  between  the  adjacent 
wells  vary  between  2.6  and  9.7  m.  Their  total  depths  vary 
between  18.2  and  20  m  below  the  land  surface  which 
occurs  between  849.32  m  and  849.64  m  above  mean  sea 
level  (AMSL)  in  the  central  area. 

[9]  Some  key  structural  information  about  the  hydrogeo¬ 
logical  units  and  properties  at  this  site  has  been  obtained 
from  neutron  (porosity)  logs  and  cross-hole  ground-pene¬ 
trating  radar  (GPR)  data.  The  neutron-log  data  were  meas¬ 
ured  every  0.06  m  in  each  borehole  and  6  values  were 
obtained  from  the  measured  count  rate  through  a  petro¬ 
physical  transform  [ Barr  ash  and  Clemo,  2002].  At  the 
BHRS,  neutron-log  data  are  primarily  sensitive  to  total  po¬ 
rosity  in  the  water-saturated,  unconsolidated,  coarse  (high- 
energy)  sedimentary  deposits  that  have  been  documented 
to  have  virtually  no  silt  and  clay  [Barrash  and  Reboulet, 
2004].  That  is,  the  measured  porosity  values  can  be  taken 
as  equivalent  to  effective  porosity  values  for  this  aquifer 
with  the  risk  of  only  very  limited  overestimation.  On  the 


basis  of  these  </>  logs,  Barrash  and  Clemo  [2002]  identified 
five  (j)  units  with  distinct  spatial  distribution,  mean,  and 
variance  (Figure  2).  The  fitted  Gaussian  functions  to  the  (f> 
distributions  of  units  1-5  have  means  equal  to  0.18,  0.24, 
0.172,  0.224,  and  0.425,  respectively,  and  standard  devia¬ 
tions  equal  to  0.022,  0.038,  0.024,  0.05,  and  0.055,  respec¬ 
tively.  Unit  5  is  a  high  cf>  channel  that  thickens  toward  the 
Boise  River  and  pinches  out  in  the  center  of  the  well  field. 
Units  1-4  are  a  sequence  of  conglomerates  with  gravel  and 
cobble  framework  and  sand  to  pebble  matrix  in  the  inter¬ 
stices  of  the  framework.  More  recently,  electrical  capacitive 
conductivity  measurements  have  allowed  the  identification 
of  a  subunit  2b,  which  is  present  in  all  of  the  measured  A-C 
wells  shown  in  Figure  1  except  wells  Bl,  B3,  Cl,  and  C2 
[Mwenifumbo  et  al.,  2009].  The  contacts  between  units  1, 
2a-2b,  3,  4,  and  5  are  projected  between  wells  B6,  Al,  and 
B3  in  Figure  2. 

[10]  Numerous  cross-hole  GPR  measurements  have  also 
been  acquired  at  the  BHRS,  and  from  these,  3 1  intersecting 
cross-hole  GPR  profiles  have  been  inverted  together  to 
obtain  a  highly  resolved  and  internally  consistent  model  of 
radar  velocity  distribution  along  the  various  profile  direc¬ 
tions  [Dafflon  et  al.,  2011].  The  obtained  velocity  distribu¬ 
tion  generally  correlates  well  with  complementary  (j)  log 
data,  with  correlation  coefficients  varying  between  —0.32 
and  —0.79  over  the  site  and  with  a  mean  equal  to  —0.57. 
Moreover,  such  GPR  velocity  tomograms  are  an  important 
source  of  </»  structural  information  between  the  boreholes 


B6  Al  B3 


Figure  2.  </>  logs  in  wells  B6,  Al,  and  B3.  For  considera¬ 

tion  of  the  general  hydrogeological  structure  at  the  BHRS, 
contacts  between  units  1-5  are  indicated  based  on  informa¬ 
tion  contained  in  the  logs  (shown  here).  Relative  horizontal 
positions  of  well  logs  are  not  to  scale. 
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because  of  their  continuous  spatial  coverage  and  the  strong 
petrophysical  relationship  (negative  correlation)  between 
GPR  velocity  and  0  in  the  aquifer  at  the  BHRS  (except  for 
unit  2b).  Figure  3a  shows  the  GPR  velocity  tomogram  along 
the  profile  B6-B3  from  the  inversion  of  3 1  profiles  together. 
The  imaged  structures  correlate  well  with  the  <j>  log  data 
(Figure  2).  Also,  in  Figure  3a,  units  1  and  3  show  relatively 
low  internal  variability,  whereas  units  2a  and  4  show  higher 
variability  and  the  presence  of  smaller  scale  structures 
(lenses).  This  is  in  agreement  with  the  previous  geostatisti- 
cal  work  of  Barrash  and  Clemo  [2002],  which  found  that 
the  variance  in  4>  in  units  2  and  4  (see  values  above)  is 
~  1.5-2  times  larger  than  that  in  units  1  and  3.  We  also  note 
that  the  GPR  clearly  imaged  the  boundary  between  units  2a 
and  2b  and  thus  corroborates  the  result  obtained  from  elec¬ 
trical  capacitive  conductivity  measurements  that  indicate 
unit  2  can  be  subdivided  into  units  with  different  ^-electric/ 
dielectric  petrophysics  [ Mwenifumbo  et  al.,  2009]. 

[11]  Finally,  3-D  realizations  of  the  <0  distribution  at  the 
BFIRS  have  been  generated  from  the  information  contained 
in  the  multidirectional  GPR  velocity  model  and  the  <j> 
logs  (B.  Dafflon  and  W.  Barrash,  3-D  stochastic  estimation 
of  porosity  distribution:  Benefits  of  using  GPR  velocity 
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Figure  3.  (a)  Velocity  tomogram  for  profile  B3-B6  that 

was  obtained  by  inverting  together  3 1  cross-hole  GPR  data 
sets  [Dafflon  et  al.,  2011].  Intersections  with  other  profiles 
are  delineated  (gray),  (b)  <0  distribution  along  profile  B3-B6, 
extracted  from  a  3-D  realization  obtained  from  the  <f>  logs 
and  velocity  tomograms  using  a  simulated-annealing  condi¬ 
tional  simulation  approach,  shown  with  the  main  lithological 
unit  boundaries  identified  at  the  BHRS.  This  realization  can 
be  considered  as  representative  of  a  large  number  of  realiza¬ 
tions  with  regard  to  implications  for  this  study. 


tomograms  in  simulated-annealing-based  or  Bayesian  se¬ 
quential  simulation  approaches,  submitted  to  Water 
Resources  Research,  2011)  through  a  simulated  annealing 
approach  [Dafflon  et  al. ,  2009] .  In  short,  this  method  allows 
the  successful  fusion  of  the  relatively  large-scale  informa¬ 
tion  contained  in  the  GPR  velocity  model  and  the  smaller- 
scale  information  contained  in  the  cf>  logs.  To  do  this,  a 
Monte  Carlo-based  process  is  used  to  optimize  the  struc¬ 
tural  variability  of  a  realization  (i.e.,  given  a  geostatistical 
target  function)  by  sequentially  and  randomly  perturbing 
this  realization,  where  each  cell  is  conditioned  on  the  infor¬ 
mation  available  at  that  cell.  This  means  that,  in  addition  to 
the  available  GPR  velocity  and  <0  data,  this  approach  uses 
the  conditional  relation  inferred  between  these  properties 
along  the  wells  and  with  geostatistical  functions  obtained 
from  the  data.  This  approach  has  been  shown  to  provide 
realistic  3-D  simulations  and  to  increase  the  accuracy  of 
cj)  prediction  compared  to  not  using  the  GPR  data  (Dafflon 
and  Barrash,  submitted  manuscript,  2011).  Thus,  the  best 
3-D  0  estimates  now  available  at  the  BHRS  are  obtained 
through  the  coupling  of  the  <0  logs  and  cross-hole  GPR  ve¬ 
locity  data.  One  representative  realization  of  the  <0  structure 
is  shown  in  Figure  3b  along  the  profile  B6-B3.  It  is  evident 
by  inspection  that  the  spatial  variability  and  the  (f>  values  in 
each  unit  are  very  consistent  with  the  </>  data  (Figure  2)  and 
the  GPR  velocity  imaged  structures  (Figure  3  a). 

3.  Tracer  Test 

[12]  The  tracer  test  considered  in  this  work  was  conducted 
at  the  BHRS  in  August  2001  [Barrash  et  al.,  2002],  A  con¬ 
servative  tracer  was  injected  into  a  4-m  interval  that  spanned 
the  boundary  between  units  2a  and  3  in  well  B3,  and  was 
removed  by  pumping  at  well  B6  from  a  4-m  interval  located 
at  the  same  elevation  as  the  injection  interval  (Figure  4). 
Wells  B3  and  B6  are  6.9  m  apart  and  their  alignment  is  ori¬ 
ented  ~30°  west  of  the  natural  head  gradient  direction.  A 
third  well  (Al)  is  located  between,  and  on  the  same  align¬ 
ment  as  wells  B3  and  B6  (Figure  1).  Well  Al  was  equipped 
to  collect  samples  from  20  0.25  m-long  sampling  zones  iso¬ 
lated  by  packers  and  situated  in  a  5-m  interval  above 
836.575  m  elevation  (Figure  4),  and  accessed  by  ports  and 
dedicated  tubing.  The  same  protocol  was  used  at  wells  Bl, 
B2,  B4,  and  B5  to  collect  samples  from  six  1-m-long  zones 
centered  on  the  4-m  interval  of  the  injection  zone,  although 
these  wells  were  on  the  margin  of  the  expected  plume  path. 
Similarly,  samples  were  taken  from  six  1-m-long  intervals  at 
B6  (pumping  well),  which  included  four  1-m-long  intervals 
that  had  perforations  in  the  riser  to  allow  pumping  intake 
over  the  injection  zone’s  elevation  span.  Flowrate  measure¬ 
ments  were  taken  with  a  digital  flowmeter  for  injection  and 
pumping  in  wells  B3  and  B6,  respectively.  (See  Barrash 
et  al.  [2002]  and  Hausrath  et  al.  [2002]  for  additional  details 
on  well  instrumentation,  sampling  protocols,  measurements, 
and  results  for  the  tracer  test.) 

[13]  The  experiment  started  with  the  injection  of  the  bro¬ 
mide  tracer  with  a  concentration  of  7.56  g  L-1  in  well  B3 
for  33.3  min  at  a  rate  of  111.6  L  min-1.  Approximately  35 
min  after  the  injection  was  completed,  pumping  was  started 
from  well  B6  with  a  low  rate  of  ~20  L  min- 1 .  These  tracer 
test  operational  parameters  were  designed  based  on  prior 
numerical  modeling  [Barrash  et  al.,  2002].  The  use  of  such 
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Figure  4.  Geometry  of  the  tracer  test  experiment  per¬ 
formed  in  2001  at  the  BHRS.  Wells  extend  above  to  land 
surface  at  ~849.5  m  elevation,  and  below  to  ~83  1 .2  m 
elevation. 


a  pumping  rate  was  intended  to  (1)  assure  the  plume  pas¬ 
sage  through  Al,  (2)  ensure  that  the  plume  would  remain 
within  the  ring  of  B  wells,  (3)  completely  remove  the  tracer 
from  the  field,  and  (4)  minimize  the  influence  of  pumping 
on  the  interval  between  wells  B3  and  Al.  Although  it  was 
planned  to  pump  continuously  at  the  same  rate  for  the 
whole  duration  of  the  test,  some  variations  occurred  due  to 
two  power  failures,  cessation  of  pumping  for  GPR  tomog¬ 
raphy,  and  some  minor  drift  and  fluctuation  of  the  pumping 
rate  [see  Barrash  et  al.,  2002,  Figure  34].  The  main  fea¬ 
tures  of  the  pumping  rate  measured  through  the  first  14.1  d 
are  as  follows:  a  rate  equal  to  ~  1 8.69  L  min~'  between 
35  min  and  6.09  ds,  then  17.66  L  min  1  until  day  14.1  with 
three  interruptions:  between  8.93  and  9.33  d,  11.81  and 
12.09  d,  and  13.81  and  14.1  d.  Finally,  after  14.1  d  the 
pumping  rate  was  increased  to  98.42  L  min-1  to  ensure  a 
nearly  complete  recovery  of  tracer  from  the  aquifer.  In  this 
regard,  >95%  of  bromide  mass  was  recovered  [Hausrath 
et  al.,  2002], 

[14]  Groundwater  samples  of  fluid  for  concentration 
measurements  were  collected  with  peristaltic  pumps  from 
the  sampling  zones  in  well  Al  prior  to  the  injection  to  get 
background  concentration  values,  from  the  mixing  tank 
prior  to  injection,  from  the  injection  line  at  4  min  intervals 
during  injection,  from  the  sampling  zones  in  Al  and  all 
B-wells  except  B3,  and  from  the  pumping  line  at  B6  every 
4  h  after  the  injection,  from  1  August  until  16  August  (see 
Hausrath  et  al.  [2002]  for  additional  detail  on  sampling 
and  measurement).  Analyses  (temperature  corrected)  with 
an  Accumet  AR50  m  and  conductivity  probe  were  com¬ 
pleted  at  the  test  site  within  4  hr  of  collection  for  all  sample 
events  during  the  test.  The  conductivity  of  the  samples  was 
measured  and  converted  to  bromide  concentration  with  a 


linear  relationship  based  on  laboratory  experiments  [ Barrash 
et  al.,  2002;  Hausrath  et  al.,  2002],  Subsequent  independent 
laboratory  tests  at  Boise  State  University  and  at  a  commer¬ 
cial  laboratory  demonstrated  that  conductivity  measurements 
and  bromide  concentrations  have  a  correlation  coefficient  of 
0.994,  and  thus  the  conductivity  measurements  and  the  cali¬ 
brated  linear  relationship  provide  close  approximation  to  the 
bromide  concentrations  over  the  full  range  of  concentrations 
observed  during  the  tracer  test  [ Hausrath  et  al.,  2002], 
Breakthrough  behavior  at  the  20  zones  in  Al  (Figure  5)  indi¬ 
cates  that  these  field  and  laboratory  methods  were  able  to 
detect  small  differences  in  tracer  concentration  between  ad¬ 
jacent  0.25-m  thick  zones  throughout  the  test.  From  Figure  5 
it  can  be  seen  that  the  concentration  breakthrough  curves  in 
well  Al  can  be  classified  into  four  vertically  contiguous 
depth  intervals  based  on  similarities  in  the  time  of  break¬ 
through  and  the  overall  magnitude  and  variation  in  the  con¬ 
centration  measurements. 

[15]  Water  level  changes  during  the  tracer  experiment 
were  mostly  related  to  important  changes  in  the  injection 
and  pumping  rates  at  the  site,  which  have  been  included 
here  in  the  modeling.  Furthermore,  the  flowrate  and  stage 
in  the  Boise  River  near  the  BHRS  were  relatively  steady. 
Daily  evapotranspiration  cycles  caused  visible  water  level 
fluctuations  in  all  of  the  surveyed  wells  with  amplitudes 
ranging  from  0.009  to  0.027  m.  These  influences  on  the 
flow  and  water  level  have  not  been  considered  in  the  hydro- 
logical  model  as  they  have  been  shown  to  not  significantly 
influence  the  results  [ Nelson ,  2007]  and  would  unnecessarily 
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Figure  5.  Tracer  test  breakthrough  curves  (BTCs)  meas¬ 
ured  in  20  sampling  intervals  along  well  Al  (shown  in  Fig¬ 
ure  4).  Based  on  their  general  behaviors  and  shapes,  the 
BTCs  have  been  divided  into  four  groups  related  to  depth 
intervals  where  measurements  were  taken. 
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complicate  the  setting  of  the  boundary  conditions  of  the 
model.  Head  change  measurements  during  the  test  were 
taken  with  vented  strain-gage  transducers  (Keller  model 
730)  in  C  and  X  wells,  prototype  fiber-optic  transducers 
(Roctest  MEMS-based  Fabray-Perot  type)  in  A  and  B 
wells,  and  e-tape  measurements  by  hand  in  X  wells  (see 
Barrash  et  al.  [2002]  for  additional  details  on  head  meas¬ 
urements  during  the  tracer  test).  Drift  after  early  time 
affected  the  prototype  fiber  optic  transducers  limiting  quan¬ 
titative  measurements  at  the  A  and  B  wells.  For  this  study 
we  decided  to  use  the  manually  measured  water  levels  in 
X-wells  only  because  they  are  the  most  accurate  measure¬ 
ments  (<0.003  m  repeatability  for  absolute  position),  and 
they  occur  at  our  model  boundary  without  imposing  much 
influence  on  local  hydrological  or  tracer  plume  behavior.  In 
fact,  these  head  measurements  are  fairly  constant  during 
each  period  of  constant  pumping  rate  (with  a  variation  in 
the  range  of  the  daily  evapotranspiration  cycles  described 
above)  and  show  some  more  clear  changes  when  the  pump¬ 
ing  rate  changes  significantly.  In  the  modeling  step,  con¬ 
stant  values  have  been  set  for  each  of  these  periods. 

[  16]  Finally,  over  the  course  of  the  tracer  test  cross-hole 
GPR  data  were  collected  almost  daily  between  B1-B4  and 
B2-B4,  and  on  three  separate  days  between  B3-B6.  For  the 
purpose  of  this  study  we  only  considered  the  data  along 
profile  B1-B4  (transverse  to  the  plume  and  intersecting 
well  Al  with  high-resolution  tracer  sampling  control),  that 
have  been  successfully  processed  to  date  [Johnson  et  al., 
2007]  and  constitute  the  most  complete  time-lapse  survey. 
To  obtain  this  profile,  data  were  acquired  with  a  100  MHz 
Mala  Geosciences  RAMAC  GPR  system  at  20  cm  spacing 
in  the  receiver  well  B 1  and  5  cm  spacing  in  the  transmitter 
well  B4.  Johnson  et  al.  [2007]  used  these  data  with  a  Fresnel 
zone  attenuations-difference  tomography  approach  (FADT) 
to  reconstruct  bulk  electrical  conductivity  changes  in  the 
subsurface  over  the  course  of  the  tracer  test.  Here  Figure  6 
shows  the  estimated  bulk  electrical  conductivity  difference 
distribution  at  various  times  along  profile  B1-B4  obtained 
by  Johnson  et  al.  [2007].  These  time-lapse  images  are  valu¬ 
able  for  semiquantitative  interpretation  of  the  large-scale 
behavior  of  the  tracer,  the  approximate  depth  reached  by 
the  sinking  tracer,  and  the  relative  lateral  spreading  or  posi¬ 
tioning  of  the  tracer. 


4.  Flow  and  Solute  Transport  Model  Setup 

[17]  Based  on  results  of  initial  modeling  of  the  BHRS 
tracer  test  [ Nelson ,  2007],  we  decided  to  model  this  test 
using  the  following  suite  of  model  programs:  MODFLOW 
[. Harbaugh  et  al.,  2000],  MT3DMS  [Zheng  and  Wang, 
1999;  Zheng,  2005],  and  SEAWAT  [e.g.,  Guo  and  Langevin, 
2002;  Langevin  et  al.,  2003,  2007],  which  are  fully  public- 
domain  codes  for  3-D  flow  and  variable-density  transport 
modeling.  SEAWAT  combines  a  modified  version  of 
MODFLOW  and  MT3DMS  into  a  single  program  that  sol¬ 
ves  the  coupled  flow  and  solute-transport  equations  to  han¬ 
dle  variable-density  transport.  This  is  done  by  iteratively 
solving  the  flow  and  transport  governing  equations  and  gen¬ 
erating  newly  calculated  heads  based  on  the  concept  of  an 
equivalent  freshwater  head  in  a  saline  groundwater  environ¬ 
ment.  That  is,  at  each  step  the  density-dependent  heads  are 
transformed  into  freshwater  heads  using  concentrations  cal¬ 
culated  in  MT3DMS  prior  to  solving  the  flow  and  transport 
equations  again. 

[is]  It  is  important  to  note  that  accounting  for  the  density 
effect  is  necessary  for  the  2001  BHRS  tracer  test  [Nelson, 
2007].  This  is  especially  evident  when  looking  at  the  shape 
and  variations  of  the  breakthrough  curves  (BTCs)  at  the 
different  sampling  intervals  in  Al  (Figure  5).  Such  sinking 
effects  have  also  been  observed  in  a  number  of  similar  field 
tracer  experiments  through  well  data  and/or  geophysical 
imaging  [e.g.,  Mackay  et  al.,  1986;  Sudicfy,  1986;  Leblanc 
et  al.,  1991 ;  Vereecken  et  al.,  2000;  Beinhorn  et  al,  2005 ; 
Muller  et  al.,  2010]  and  are  discussed  in  review  publications 
about  variable-density  flow  issues  [e.g.,  Simmons  et  al., 
2001;  Diersch  and  Kolditz,  2002;  Beinhorn  et  al.,  2005; 
Simmons,  2005],  In  the  present  study,  the  bromide  concen¬ 
tration  was  ~7.5  g  L_1  at  the  time  of  injection,  or  about  one 
fifth  of  the  concentration  of  salt  in  ocean  water,  which  was 
shown  to  be  very  significant  with  regard  to  density  effects 
[Nelson,  2007].  Enumeration  of  influences  of  tracer  concen¬ 
tration  on  sinking  and  of  concentrations  used  in  various  field 
experiments  can  be  found  in  the  work  of  Beinhorn  et  al. 
[2005].  Unfortunately,  the  inclusion  of  density  effects  in 
modeling  such  tracer  experiments  results  in  a  significant 
increase  in  computing  time  of  at  least  one  order  of  magni¬ 
tude.  Although  the  transport  can  be  simulated  in  a  subdo¬ 
main  of  the  flow  model  (i.e.,  by  identifying  specified  heads 


Figure  6.  Bulk  electrical  conductivity  difference  along  profile  B1-B4  obtained  from  time-lapse  GPR 
cross-hole  measurements  using  a  Fresnel  zone  attenuation-difference  tomography  approach  [Johnson 
et  al.,  2007]. 
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along  marginal  flow  lines)  there  is  still  a  trade-off  between 
increased  model  complexity  (because  of  nonlinear  density- 
dependent  flow)  and  computing  time. 

[19]  The  model  and  boundary  conditions  were  defined 
using  prior  information  for  this  site  and  details  of  the  tracer 
test.  Boundaries  for  the  unconfined  aquifer  model  follow 
the  X  ring  of  wells  (Figure  1)  and  thus  are  controlled  by 
head  measurements  at  the  X-wells  during  the  experiment. 
Given  the  relatively  large  distance  between  the  B-wells  and 
the  X-wells,  and  the  relatively  small  (observed)  head  varia¬ 
tion  in  the  X-wells  (section  3),  the  specified  heads  at  the 
X-wells  do  not  significantly  affect  the  flow  and  solute 
transport  behavior  locally  in  the  region  of  tracer  transport 
[Nelson,  2007],  The  natural  head  gradient  before  the  start 
of  the  experiment  was  ~0.001  in  the  flow  direction  oriented 
~30°  east  of  the  B3-B6  alignment.  Within  the  numerical 
model,  cell  dimensions  are  generally  0.2  m  x  0.2  m  x 
0.2  m  in  the  region  of  the  B-wells  (where  the  tracer  trans¬ 
port  occurs),  and  then  progressively  expand  outward  with  a 
maximum  size  oflmx  lmx0.5m.  Relevant  injection 
and  pumping  variations  during  the  test  are  included  in  the 
model.  The  injection  and  pumping  intervals  were  modeled 
with  high  K  equal  to  1  m  s~ 1  and  4>  equal  to  1  to  simulate 
well  hydraulics.  Aquifer  specific  yield,  specific  storage, 
and  the  diffusion  coefficient,  which  are  relatively  insensi¬ 
tive  parameters  in  this  experiment,  are  set  to  constants: 
Sy  =  0.21,  Ss  =  4.5e-5  m_1,  and  D  =  1.34e-9  m2  s-1.  Kis 
treated  as  isotropic  on  the  basis  of  prior  information  [Bar- 
rash  et  al,  2006].  Also,  one  parameter  that  is  included  in 
SEAWAT  and  which  influences  density  effects  is  the  slope 
of  the  linear  equation  of  state  that  relates  fluid  density  to 
solute  concentration  [Langevin  et  al.,  2007].  This  parame¬ 
ter,  commonly  set  to  0.7143  for  fresh-  and  saltwater  inter¬ 
actions,  is  approximated  by  dividing  the  density  difference 
over  the  range  of  end-member  fluids  by  the  difference  in 
concentration  between  the  end-member  fluids  [Langevin 
et  al.,  2007].  This  means  that  this  parameter  depends,  in 
part,  on  chemical  properties,  pressure,  and  temperature  of 
the  injected  tracer  solution  and  the  ambient  water.  Although 
this  slope  value  is  not  discussed  much  in  the  literature  and 
is  relatively  complex  to  estimate  in  a  field  experiment,  pre¬ 
liminary  tests  have  shown  that  changes  in  the  result  of  this 
study  are  limited  with  regard  to  the  uncertainty  involved  in 
this  parameter.  Numerical  and  laboratory  study  about  this 
parameter  would  definitely  be  worthwhile. 

[20]  Finally,  longitudinal,  horizontal  transverse,  and  verti¬ 
cal  transverse  dispersivities  have  been  set  to  0.06  m  (+0.04, 
-0.02),  0.18  m  (+0.07,  -0.08),  and  0.02  m  (+0.03, 
—0.01),  respectively,  based  on  (1)  grid  search  evaluations 
using  a  homogeneous  model,  (2)  consideration  of  the  shape 
of  the  simulated  breakthrough  curves  (BTCs)  compared  to 
the  measured  BTCs,  and  (3)  plume  position  information  in 
GPR  time-lapse  images.  The  values  in  parentheses  indicate 
ranges  of  dispersivity  values  that  influence  the  simulated 
concentration  only  by  slightly  shifting  absolute  concentra¬ 
tion  values  by  about  a  proportional  amount.  Such  variations 
have  been  shown  to  be  similar  to  those  that  could  be  pro¬ 
duced  by  error  in  positions  within  wells  (i.e.,  deviation). 
One  should  note  that,  based  on  some  additional  tests  and  on 
the  results  of  a  study  that  included  evaluating  the  effect  of 
longitudinal  and  transverse  dispersivity  values  on  the  shape 
of  a  plume  [Sudicky  et  al.,  2010],  the  influence  of  varying 


dispersivity  values  to  some  degree  (as  above),  does  not 
significantly  influence  the  shape  of  the  plume.  Although  in 
this  study  the  longitudinal  dispersivity  is  smaller  than  the 
transverse  dispersivity,  which  is  not  commonly  observed, 
we  kept  these  values  because  that  dispersivity  is  not  only 
simulating  the  spreading  of  the  plume  produced  by  nonmod- 
eled  heterogeneity  but  it  is  also  a  means  of  reproducing 
some  effects  of  parameters  that  we  did  not  include  or  failed 
to  reproduce  in  the  hydrological  model  [see  also  Sudicky 
et  al.,  2010],  such  as  (1)  spatial  trends  of  relatively  lower 
heads  further  from  the  river  associated  with  small  river  stage 
increases  and/or  evapotranspiration  [Malania  and  Johnson, 
2010;  Johnson,  2011]  that  could  produce  transverse  move¬ 
ment  away  from  the  river  and  thus  perpendicular  to  the 
main  tracer  path  direction,  (2)  possible  northeastward-slop¬ 
ing  topography  on  the  unit  2-unit  3  contact  in  the  northeast¬ 
ern  portion  of  the  B-ring  of  wells  [Barrash  et  al.,  2002],  and 
(3)  a  limitation  in  accuracy  of  modeling  the  hydrological 
processes  and  parameters  of  influence  such  as,  for  example, 
the  coupling  between  spreading  near  injection  zone  and 
near-well  hydrogeologic  structures. 

5.  Parameter  Estimation  Strategies 

[21]  We  estimated  the  distribution  of  K  and  (l>  from  tracer 
BTCs  using  three  different  strategies.  This  allows  for  com¬ 
parison  of  both  the  results  and  the  methods  for  obtaining 
them.  The  variability  of  the  obtained  results  may  influence 
our  confidence  in  the  estimated  distribution  of  hydrological 
properties.  Another  reason  to  do  this  is  to  evaluate  the 
advantages  and  drawbacks  of  various  approaches,  which 
may  differ  given  the  objective  of  a  study,  the  number  of  pa¬ 
rameters  to  estimate,  available  computing  facilities,  and 
details  of  the  forward  model.  Below,  to  evaluate  model  mis¬ 
fit  we  used  the  root-mean-square  (RMS)  difference  between 
simulated  and  measured  BTCs  using  all  available  points  in 
time,  which  have  a  relatively  constant  sampling  time  inter¬ 
val.  Although  we  tested  other  misfit  evaluation  techniques 
(first  and  second  moments,  quartiles,  and  first  arrival  time 
with  areas  below  the  BTCs),  they  did  not  significantly 
change  the  overall  results. 

5.1.  Grid  Search  Strategy 

[22]  We  first  used  a  grid  search  strategy  to  evaluate  the 
main  properties  involved  in  the  transport  process  and  to  per¬ 
form  sensitivity  analysis.  The  grid  search  approach  is  based 
on  the  evaluation  of  one  or  more  parameters  through  a  range 
of  values,  with  the  possibility  of  starting  with  a  large-scale 
search  grid  and  subsequently  refining  it.  This  means  the 
selected  search  domain,  which  is  defined  by  a  range  of  val¬ 
ues  that  we  assign  to  each  parameter,  is  sampled  fully  for 
each  dimension  of  the  problem  without  optimization,  so  this 
is  the  most  intensive  computing  method. 

[23]  The  main  advantages  of  the  grid  search  approach 
are  that  (1)  the  process  is  easily  parallelized,  even  on  a  net¬ 
work  of  independent  computers  (or  distributed  memory 
cluster)  as  no  optimization  process  is  performed;  (2)  the 
objective  function  can  be  evaluated  after  the  simulation 
process  has  been  completed,  which  in  turn,  allows  for  eval¬ 
uation  of  multiple  objective  functions;  (3)  this  approach 
does  not  require  optimization,  regularization,  or  lineariza¬ 
tion,  and  thus  is  not  sensitive  to  convergence  problems  that 
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may  be  significant  in  large,  complex  models;  and  (4)  the 
sensitivity  of  each  parameter  can  be  evaluated  because  the 
domain  of  the  solution  is  directly  sampled.  However,  a  sim¬ 
ple  visual  observation  of  the  solution  domain  is  only  possi¬ 
ble  with  a  very  small  number  of  parameters.  Finally,  we 
note  that  the  computing  efficiency  of  this  approach  can  be 
improved  using  a  Monte  Carlo  strategy  to  avoid  visiting  all 
the  cells  of  the  grid,  or  a  multistep  approach  where  only 
a  small  number  of  parameters  are  evaluated  sequentially 
and/or  the  grid  is  iteratively  redefined  [e.g.,  Stein,  1987; 
Sambridge  and Mosegaard,  2002], 

5.2.  Simulation  of  Hydrological  Property  Distributions 

[24]  Stochastic  simulations  are  used  for  a  number  of  pur¬ 
poses  in  earth  sciences  including  heterogeneity  reconstruc¬ 
tion,  uncertainty  assessment,  and  optimization  [e.g.,  Deutsch 
and  Journel,  1998;  de  Marsily  et  al.,  2005].  Such  simula¬ 
tions  can  be  unconditional  or  conditioned  to  available  data  at 
some  locations.  Although  numerous  approaches  exist  [e.g., 
Deutsch  and  Journel,  1998;  Christakos,  2000;  Tarantola, 
2005],  they  are  similar  in  that  they  require  some  statistical 
information  like  correlation  functions,  mean  and  variance,  or 
some  data  to  constrain  these  values.  Here  we  use  the  well- 
known  sequential  Gaussian  simulation  (SGS)  algorithm 
from  the  Geostatistical  Software  Library  (GSLIB)  [. Deutsch 
and  Journel,  1998],  and  we  use  information  available  at  the 
BHRS  to  define  statistical  parameters  in  the  simulation  pro¬ 
cess.  We  use  SGS  to  generate  a  large  number  of  uncondi¬ 
tional  simulations  of  hydrological  property  distributions 
(K  and  <j>,  which  are  assumed  to  be  positively  correlated  and 
linearly  related)  and  to  investigate  the  commonalities  among 
those  realizations,  which  reproduce  the  measured  BTCs  rea¬ 
sonably  well.  For  clarity,  we  note  that  simulations  used  here 
to  find  models  that  fit  the  BTCs  are  not  conditioned  to  any 
data  point,  whereas  simulations  conditioned  to  measured 
hydrological  data  are  more  commonly  performed  to  assess 
the  uncertainty  in  predicting  hydrological  processes  such  as 
flow  and/or  transport. 

[25]  In  particular,  this  procedure  generates  and  supports 
investigation  of  heterogeneity  at  a  smaller  scale  than  achiev¬ 
able  with  other  parameter  estimation  approaches.  Genera¬ 
tion  of  small-scale  heterogeneity  of  hydrological  properties 
may  support  (1)  improving  fits  to  BTCs  and  (2)  finding 
models  showing  the  nonuniqueness  in  the  solution  domain 
(i.e.,  finding  models  that  fit  the  data  to  a  similar  level  and 
evaluating  how  different  they  are).  The  only  critical  issue 
related  to  this  approach  is  the  need  to  reduce  computing 
time  by  constraining  the  simulation  process  with  prior  sta¬ 
tistical  information,  such  as  the  global  mean  and  variance  of 
the  property  distributions,  which  limit  the  search  domain 
and  thus  may  influence  the  results. 

5.3.  Regularized  Least  Square  Inversion 

[26]  Approaches  for  the  inversion  of  hydrological  flow 
and/or  solute  transport  consist  of  sequentially  updating  the 
parameters  in  a  model  to  minimize  the  difference  between 
the  calculated  and  observed  data.  Reviews  of  the  state  of  the 
art  of  such  techniques  are  numerous  [e.g.,  Menke,  1984; 
Yeh,  1986;  de  Marsily  et  al.,  1999;  Aster  et  al.,  2005; 
Carrera  et  al.,  2005;  Tarantola,  2005],  Here  we  apply  a 
deterministic  inversion  that  optimizes  the  value  of  parame¬ 
ters  at  pilot  points  [e.g.,  Doherty,  2003;  Carrera  et  al.. 


2005;  Johnson  et  al.,  2009]  subject  to  an  objective  function 
that  contains  a  least  squares  data  misfit  criterion  and  a  regu¬ 
larization  term.  The  objective  function  is 

$Kst)  =  \\Wc(Ge[mest]  -  Cobs)||2  +  0\\Wmmesl\\2 ,  G) 

where  Gc(mest )  and  C0bs  are  vectors  of  length  n,  where  n  is 
the  number  of  measurements.  Gc(mest )  contains  the  concen¬ 
tration  data  predicted  from  modeling  the  solute  transport 
(using  SEAWAT)  given  the  estimated  parameters  mest  and 
Cobs  as  the  measured  concentration  data.  Wc  is  the  data 
weighting  matrix  of  size  n  x  n,  where  each  diagonal  ele¬ 
ment  is  the  reciprocal  of  the  estimated  standard  deviation 
of  each  measurement  and  the  off-diagonal  elements  are 
zero  if  data  errors  are  uncorrelated.  The  second  term  of 
equation  (1)  is  the  regularization  term,  where  Wm  (of  size 
m  x  m,  where  m  is  the  number  of  unknowns)  allows  direc¬ 
tional  smoothness  constraints  and  0  is  a  scalar  weighting 
parameter  that  defines  how  much  the  regularization  influen¬ 
ces  the  process  of  decreasing  the  misfit  between  the  calcu¬ 
lated  and  measured  data;  Wm  and  0  can  be  directionally 
dependent. 

[27]  At  each  iteration,  we  update  parameters  by  finding  a 
step  A»2est  that  reduces  the  value  of  the  objective  function. 
After  developing  equation  (1)  through  a  Taylor  series 
approximation,  the  final  equation  is 

(JTCWTCWCJC  +  0WTm  Wm)Amest 

(2) 

=  JTo  W0Wc(Cobs  -  Gc [»7?est] )  -  0WTm  Wmmest, 

where  the  concentration  Jacobian  matrix  Jc  is  obtained 
using  a  finite  difference  approach  [Johnson  et  al.,  2009]. 
The  perturbation  amount  A mest  is  solved  for  each  iteration 
using  a  conjugate  gradient  least  squares  algorithm.  After 
each  iteration,  the  new  estimated  mest  is  obtained  by  adding 
to  the  previous  solution  the  updated  vector  Amest,  which  is 
rescaled  on  the  basis  of  finding  the  optimum  scaling  that 
produces  the  greatest  decrease  to  the  objective  function. 
The  ideal  number  of  processors  to  solve  such  a  problem  is 
one  plus  the  number  of  pilot  points  for  which  the  sensitivity 
matrix  needs  to  be  evaluated.  In  this  study  we  used  an 
eight-core  cluster  and  the  parallel  computation  tools  avail¬ 
able  in  MATLAB. 

[28]  Also,  in  this  study,  as  in  hydrological  inversion  in 
general,  the  vector  mest  contains  the  estimated  hydraulic  pa¬ 
rameters  at  a  set  of  pilot  points  only.  These  values  distrib¬ 
uted  at  a  few  cells  in  the  model  domain  are  then  used  to 
generate  values  at  each  grid  element.  Generally,  this  is 
done  through  interpolation,  kriging,  and/or  simulation  proc¬ 
esses  [e.g.,  Certes  and  Demarsily,  1991;  Ramarao  et  al., 
1995;  Rubin  et  al.,  2010],  The  choice  of  the  appropriate 
method  depends  on  the  number  of  pilot  points  and  some 
prior  idea  of  the  required  heterogeneity  to  fit  the  measured 
data,  which  can  be  gained  by  preliminary  modeling  and 
data  mining,  and/or  from  structural  information  available 
from  other  geophysical  or  hydrological  data  at  the  site.  The 
more  pilot  points  that  are  used  to  characterize  a  spatial  pa¬ 
rameter  distribution,  the  less  important  is  the  choice  of 
interpolation  method  [Doherty,  2003].  Also,  it  is  important 
to  note  that  kriging  the  values  from  the  pilot  points, 
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although  producing  some  smoothness  in  the  model,  does 
not  control  the  correlation  between  the  pilot  points.  In  fact, 
the  overall  variability  of  the  parameters  defined  at  the  pilot 
points  is  controlled  by  the  regularization  term  only. 

[29]  Although  the  implementation  of  this  regularized 
inversion  approach  looks  more  complex  than  previously  dis¬ 
cussed  methods,  once  the  equation  system  is  set,  advantages 
include  (1)  the  significant  reduction  of  computing  resources 
(i.e.,  although  it  is  less  easy  to  parallelize  on  a  network  of 
independent  computers),  (2)  the  ability  to  find  a  model  of 
minimal  variability  that  reproduces  the  data,  and  (3)  the 
ability  to  solve  for  a  relatively  large  number  of  parameters 
in  general.  In  particular,  the  reduction  of  computing  time 
allows  us  to  more  easily  evaluate  various  parameterizations 
and/or  problem  configurations  and  their  effect  on  the  result¬ 
ing  distribution  of  hydrological  properties. 

6.  Results 

6.1.  Parameter  Estimation  Through  a  Grid  Search 
Strategy 

[30]  The  tracer  test  data  are  first  investigated  with  a  grid 
search  strategy  to  evaluate  some  simple  model  configurations 
(i.e.,  defined  by  a  relatively  small  number  of  parameters)  that 
may  allow  us  to  reproduce  these  data  to  a  first  order.  We  first 
evaluate  a  homogeneous  model,  which,  because  a  limited 
number  of  parameters  are  evaluated,  is  relatively  unifonnly 
sampled  and  the  residual  misfit  between  predicted  and  meas¬ 
ured  concentration  data  can  be  observed  throughout  the  full 
considered  parameter  domain.  Then  we  investigate  how  a 
relatively  simple  layered  model  can  fit  measured  data  and 
better  predict  measured  concentrations  than  a  homogeneous 
model. 

6.1.1.  Homogeneous  Model 

[31]  Here  we  evaluate  some  broad  ranges  of  <l>  and  K  for 
a  homogeneous  model.  Figure  7  shows  RMS  residuals 
between  all  predicted  and  measured  concentrations  for  each 
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Figure  7.  Grid  search  results  for  the  evaluation  of  homo¬ 
geneous  models  of  </>  and  K  that  support  the  prediction  of 
solute  concentration  data  and  are  considered  here  in  regard 
to  the  root-mean-square  (RMS)  difference  between  pre¬ 
dicted  and  measured  BTCs  at  well  Al. 


set  of  evaluated  properties.  It  is  clear  that  relatively  similar 
fits  result  from  correlated  changes  in  <j>  and  K.  This  is 
explained  by  the  link  between  4>  and  K  in  the  transport 
equations.  However,  the  RMS  residual  increases  signifi¬ 
cantly  for  cj)  lower  than  ~0.28,  so  it  appears  that  the  cou¬ 
pling  between  (j)  and  K  has  some  influence  on  the  sinking 
(density  effect)  of  the  tracer,  and  the  4>  and  K  values  need  to 
be  in  a  limited  range  to  fit  the  BTCs  acceptably.  It  is  impor¬ 
tant  to  note  that  only  enforcing  very  high  K  values  at  some 
depth  interval  does  not  allow  draining  of  the  tracer  and 
improved  fitting  to  the  BTCs.  Instead,  both  high  <f>  and  K 
seem  to  be  necessary  in  order  for  a  larger  amount/mass  of 
tracer  to  sink  rapidly  such  that  adequate  fits  to  BTCs  are 
obtained. 

[32]  Figure  8  shows  the  predicted  concentration  from  the 
(j)  and  K  model  with  the  best  fit  of  the  BTCs.  Although  this 
model  gives  acceptable  fits  to  first  arrival  time  and  some 
variations  in  the  BTCs,  and  thus  is  worthy  for  a  first-order 
evaluation  of  hydrological  data  and  properties,  more  hetero¬ 
geneity  is  needed  to  get  better  fits  at  each  measurement 
depth  level.  For  example,  the  velocity  in  the  bottom  part  of 
the  model  should  be  larger  to  better  fit  the  early  arrival  time 
between  836.7  m  and  837.45  m.  Also,  the  homogeneous 
model  does  not  allow  the  tracer  to  sink  enough  and  thus  pre¬ 
dicted  concentrations  between  838.45  m  and  839.95  m  are 
too  high  and  the  breakthrough  arrival  time  is  too  early  there. 

6.1.2.  Layered  Model 

[33]  A  grid  search  approach  is  used  to  evaluate  a  three¬ 
layered  model  with  various  possible  layer  thicknesses  and 
K  and  <p  values.  For  each  layer,  the  ranges  of  values  eval¬ 
uated  for  these  parameters  were  defined  using  prior  infor¬ 
mation  from  a  GPR  velocity  model  \Dafflon  et  al.,  2011] 
and  <j>  logs  [Barrash  and  Clemo,  2002],  and  from  the  grid 
search  using  a  homogeneous  model.  Thus,  </>,  K,  and  the  top 
of  the  layer  are  evaluated  for  the  interval  [0.15;  0.25;  0.35], 
[0.5e-4;  le-4;  3e-4],  and  [850],  respectively,  for  the  upper 
layer;  [0.25;  0.35;  0.45],  [2e-4;  4e-4;  6e-4],  and  [838.475; 
838.075;  837.675],  respectively,  for  the  central  layer;  and 
[0.25;  0.35;  0.45],  [2e-4;  4.5e-4;  7e-4],  and  [837.275; 
836.675;  836.075],  respectively,  for  the  lower  layer. 

[34]  The  reason  for  evaluating  a  three-layered  model  and 
such  a  small  number  of  parameters  and  range  of  values  is 
the  significant  demand  on  computing  resources,  which  in  this 
case  was  ~20  dual-core  computers  running  at  night  for  1 
week  to  perform  6561  simulations.  Computing  time  would 
increase  exponentially  with  the  number  of  evaluated  parame¬ 
ters,  such  as,  for  example,  evaluating  various  dispersivity  val¬ 
ues.  Hydrologically,  use  of  a  three-layered  model  rather  than 
a  five-layered  model  for  the  BHRS  (e.g.,  Barrash  and  Clemo 
[2002];  and  see  Figures  2  and  3  in  this  work)  is  justified 
because  flow  and  transport  during  the  tracer  test  occurred  pri¬ 
marily  or  exclusively  in  the  middle  layers  of  the  system. 

[35]  Figure  9  shows  the  model  with  the  lowest  RMS  dif¬ 
ference  (0.0421)  between  predicted  and  measured  concen¬ 
trations  among  all  of  the  evaluated  models.  The  fit  to  BTCs 
is  improved  compared  to  Figure  8,  mainly  because  K  and  0 
are  now  larger  below  838.075  m  (also  see  Figure  2).  In 
fact,  the  ratio  of  K  to  </>  is  higher  which  implies  higher  ve¬ 
locity  and  earlier  arrival  of  large  concentrations  below 
837.45  m.  The  arrival  times  and  principal  behaviors  of  the 
BTCs  are  relatively  well  reproduced.  The  concentrations 
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Figure  8.  Best  fit  model  obtained  through  a  grid  search  evaluation  of  homogeneous  K  and  (j>  fields : 
(a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown  in  Figure 
5),  (b)  estimated  <j>,  and  (c)  estimated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e)  simulated  concen¬ 
tration  after  9.3  and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 


simulated  between  B 1-B4  (Figure  9e)  show  the  presence  of 
the  tracer  at  depths  that  are  similar  to  those  in  time-lapse 
GPR  images  (Figure  6).  Also,  the  significant  changes  in  <f> 
and  K  at  838.075  m  are  consistent  with  the  sharp  boundary 
between  units  2a  and  2b  at  838  m  in  the  GPR  tomogram 
along  profile  B6-A1-B3  (Figure  3)  and  in  the  <f>  logs  where 
high-porosity  lenses  >0.30  occur  in,  and  likely  between, 


wells  B3  and  A1  between  838  and  836  m  elevation.  Indeed, 
the  results  indicate  that  a  simple  three-layer  model  can 
explain  a  large  part  of  the  information  contained  in  the  data, 
and  thus  that  principal  behaviors  in  the  BTCs  can  be  related 
to  hydrostratigraphic  layers  or  large  hydrological  struc¬ 
tures.  However,  some  behaviors  in  the  measured  BTCs  that 
are  still  not  well  reproduced  by  the  model,  such  as  the 
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Figure  9.  Best  fit  model  obtained  through  a  grid-search  evaluation  of  a  three-layer  model  with  move- 
able  boundaries:  (a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along 
A1  (shown  in  Figure  5),  (b)  estimated  </>,  and  (c)  estimated  K  distributions,  shown  along  profile  B6-B3 ; 
(d-e)  simulated  concentration  after  9.3  and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4, 
respectively. 
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decrease  in  concentration  after  ~10  d  in  the  BTCs  measured 
between  836.7  and  837.45  m,  will  be  investigated  next  by 
simulating  small-scale  variability  in  the  hydrological  models. 

6.2.  Stochastic  Simulation  of  Hydrological  Properties 

[36]  We  perform  stochastic  simulations  to  obtain  realiza¬ 
tions  of  K  and  linearly  related  <fi  with  small-scale  variability 
that  may  improve  matches  to  the  measured  BTCs.  We 
investigate  this  with  two  approaches:  (1)  with  the  three- 
layer  model  shown  in  Figure  9  (i.e.,  an  a  priori  setting  of  the 
large-scale  structures)  to  examine  improvement  in  fitting 
the  BTCs  due  to  the  simulation  of  small-scale  heterogeneity 
in  each  layer,  and  (2)  with  simulations  that  reproduce  small- 
and  large-scale  structures  without  prior  layering  specified  in 
the  model  to  evaluate  the  extent  to  which  such  simulations 
show  similarities  when  reproducing  the  measured  data  to 
some  degree.  It  is  important  to  note  that  parameterizations 
of  these  two  simulation  approaches,  involving  means  and 
covariance  functions,  are  different  due  to  the  different  scale 
at  which  heterogeneity  needs  to  be  defined  in  each. 

6.2.1.  Unconditional  Simulations  in  a  Three-Layer 
Model 

[37]  From  the  three-layer  models  with  the  best  fit  BTCs 
from  the  grid  search  estimation  process  (section  6.1.2),  an 
arbitrary  number  of  2500  simulations  were  performed  in 
each  layer  to  create  small-scale  heterogeneity.  The  mean 
and  variance  of  the  simulations  in  each  layer  were  defined 
based  on  the  results  obtained  earlier  through  the  grid  search 
strategy  (Figure  9).  For  the  bottom,  middle,  and  top  layers, 
mean  K  is  set  to  5.3e-4,  4.2e-4,  and  2e-4  m  s-1,  respectively 
(i.e.,  similar  to  values  from  fully  penetrating  pumping  tests 
at  the  BHRS  [Barrash  et  al.,  2006]),  with  standard  devia¬ 
tions  of  2e-4,  2e-4,  and  le-4,  respectively,  and  the  minimum 


value  is  set  to  0.5e-4  m  s_I.  Also,  <t>  is  linked  to  K  with  a 
linear  relationship,  which  is  consistent  with  the  positive  cor¬ 
relation  observed  between  </>  and  K  in  the  three-layer  model, 
and  with  the  need  to  set  some  parameters  using  prior  results 
to  simplify  the  problem.  For  the  bottom,  middle,  and  top 
layers,  the  resulting  mean  <j>  values  are  0.33,  0.33,  and  0.2, 
respectively,  with  standard  deviations  of  0.08,  0.08,  and  0.08, 
respectively,  and  the  minimum  value  is  set  to  0.05.  The  cor¬ 
relation  lengths  were  defined  based  on  the  vertical  correlation 
length  of  1 .2  m  obtained  from  fitting  the  correlation  function 
from  the  (j>  logs  for  the  various  depth  intervals,  which  in  turn, 
is  similar  to  values  determined  by  Barrash  and  Clemo 
[2002].  Other  values  could  be  used  to  parameterize  the  simu¬ 
lations,  but  these  are  reasonable  choices  and  the  purpose  here 
is  more  to  evaluate  improvement  in  the  BTCs  by  introducing 
small-scale  structures  than  to  estimate  all  possible  solutions. 

[38]  Figure  10  shows  the  hydrological  realization  with 
the  best  fit  to  the  measured  BTCs  (RMS  residual  is  0.0248). 
This  simulation  shows  the  extent  to  which  small-scale  struc¬ 
ture  can  improve  fits  to  the  BTCs,  including  reproduction  of 
some  small-scale  behaviors  in  the  BTCs.  The  five  “best” 
simulations  have  residuals  between  0.0248  and  0.0287  and 
the  RMS  residual  distribution  of  all  the  simulations  has  a 
Gaussian  behavior  that  can  be  fit  with  mean  and  standard 
deviations  of  0.047  and  0.008,  respectively.  This  shows 
that  a  relatively  small  number  of  solutions  significantly 
improve  the  three-layer  model  (Figure  9)  and  that  the  simu¬ 
lation  of  such  small-scale  structures  can  strongly  influence 
behaviors  in  the  predicted  concentration  curves  (Figure 
10).  Also,  we  note  that  high  cj)  values  were  required  to 
improve  the  fits,  which  again,  are  consistent  with  the  local 
structure  of  high  cj)  lenses  below  838  m  elevation  (i.e.,  espe¬ 
cially  between  838  and  836  m  but  also  between  836  and 


a) 

s  °-4 

u  0.2 

0 


0 


10 


b)  B6 
842 


CIC0 


15 


46  48  50  52 

X[m]  r  ., 

B6  Al  B3^[m/s] 


48  50  52  48  50  52 
.  X[m]  X[m] 

e)  B1  Al  B4  B1  Al  B4 


C/C0 


5  10 

time  [days] 

;  RMS  difference  =  0.0248 


46  48  50  52 
X[m] 


2  4 
D[m] 


Figure  10.  Best  fit  model  obtained  through  a  simulation  approach  using  the  geometry  of  the  three-layer 
model  shown  in  Figure  9:  (a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals 
along  Al  (shown  in  Figure  5),  (b)  simulated  <j>,  and  (c)  simulated  K  distributions,  shown  along  profile 
B6-B3;  (d-e)  simulated  concentration  after  9.3  and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4, 
respectively. 
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834  m;  see  Figure  2)  and  corroborates  the  results  of  the 
grid  search  approach.  Finally,  although  the  model  shown  in 
Figure  10  gives  a  very  good  fit  to  the  BTCs  overall,  some 
features  are  still  not  well  reproduced,  especially  some  high- 
frequency  variations  in  the  BTCs,  such  as,  for  example,  the 
decrease  in  concentration  after  ~10  d  and  then  the  increase 
again  after  ^12  d  in  the  BTCs  measured  between  836.7 
and  837.45  m  elevation  (Figure  10).  In  fact,  incomplete  fit¬ 
ting  of  BTCs  is  common  [e.g.,  Hyndman  etal.,  2000;  Scheibe 
and  Chien,  2003;  Johnson  et  al.,  2009]  and  it  is  difficult  to 
predict  the  degree  to  which  real  BTCs  can  be  reproduced. 
This  is  mainly  related  to  limitations  in  the  understanding  and 
modeling  of  transport  processes,  boundary  effects,  and  the  pa¬ 
rameters  of  influence. 

6.2.2.  Unconditional  Simulations  in  a  One-Layer 
Model 

[39]  Now  we  investigate  if  structural  similarities  can  be 
found  between  realizations  that  fit  the  BTCs  relatively  well 
but  have  not  been  conditioned  with  prior  layer  information. 
To  this  end,  simulations  of  K  and  linearly  related  4>  fields 
are  performed  based  on  previously  estimated  global  statis¬ 
tics  at  the  BHRS.  The  global  means  of  K  and  <l>  for  the  sim¬ 
ulations  are  set  to  4e-4  m  s'1  and  0.3,  respectively,  and 
their  standard  deviations  are  set  to  2e-4  and  0.1,  respec¬ 
tively.  These  are  similar  to  thickness-averaged  K  values 
from  fully  penetrating  pumping  tests  [Barrash  et  al.,  2006] 
and  important  <fi  structures  at  B3  and  Al  (Figure  2).  Mini¬ 
mum  (j)  and  K  values  are  set  to  0.05  and  0.5e-4,  respec¬ 
tively.  Florizontal  and  vertical  correlation  lengths  have 
been  set  to  8  m  and  1.6  m,  respectively.  These  lengths  are 
based  on  the  previous  results  including  structural  informa¬ 
tion  from  GPR  data  (Figure  3),  and  are  similar  to  lengths 


obtained  from  the  geostatistical  analysis  of  K  from  slug  test 
results  at  the  BHRS  [Cardiff et  al.,  2011]. 

[40]  The  model  with  the  lowest  RMS  difference  between 
predicted  and  measured  BTCs  is  shown  in  Figure  11,  and 
Figure  12  shows  four  other  realizations  with  next  best  fits  to 
the  BTCs.  These  five  models  have  RMS  residuals  between 
0.026  and  0.0276  and  show  structural  similarities  with  (1) 
layers  of  higher  K  and  ^  between  836  and  838  m  (Figures 
11,  12c,  and  2),  which  also  are  visible  in  the  mean  of  these 
simulations  (Figure  12b);  and  (2)  the  strong  change  to  lower 
(j>  and  K  above  838  m  (that  also  occurred  in  the  three-layer 
model;  see  Figures  9  and  10)  compared  to  the  strong  change 
in  GPR  velocity  and  </>  at  the  same  elevation  (Figures  2  and 
3a).  Compared  to  the  three-layer  model  obtained  by  grid 
search  (Figure  9),  the  BTCs  are  better  fit  with  the  lowest 
RMS  models  (Figures  11  and  12c),  which  is  due  to  more 
structural  variability  in  the  simulations.  In  particular,  the 
choice  of  a  relatively  large  lateral  correlation  length  com¬ 
pared  to  the  size  of  the  model  allows  some  variation  in  the 
mean  or  variance  locally,  which  in  turn,  allows  a  large  dis¬ 
tribution  of  RMS  differences  (Figure  12a)  and  improves  the 
possibility  of  finding  simulations  with  particular  features 
that  are  able  to  fit  the  data  relatively  well  (Figures  1 1  and 
12c),  even  if  such  simulations  are  few. 

6.3.  Parameter  Estimation  Through  a  Regularized 
Inversion 

[41]  In  section  6.3  we  perform  regularized  inversion  with  0 
as  (1)  a  predefined  constant  value  (section  6.3.1),  (2)  an  esti¬ 
mated  parameter  in  the  inversion  process  (section  6.3.2),  and 
(3)  a  predefined  distribution  obtained  from  3-D  (j>  recon¬ 
struction  from  geophysical  data  (section  6.3.3).  Such 
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Figure  11.  Best  fit  model  obtained  through  a  simulation  approach  without  layers:  (a)  measured  (thin 
lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  Al  (shown  in  Figure  5),  (b)  simulated 
<f),  and  (c)  simulated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e)  simulated  concentration  after  9.3 
and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 
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Figure  12.  (a)  Distribution  of  the  root-mean-square  difference  between  measured  and  predicted  con¬ 

centrations  obtained  from  6000  simulations  of  the  hydrological  structure,  (  b)  Mean  of  the  linearly  related 
K  and  4>  realizations  of  the  five  best  fits  to  the  measured  BTCs  (red  circle  in  Figure  12a).  The  best-fitting 
realization  is  shown  in  Figure  11.  (c)  The  next  four  best  fits,  (d)  Four  realizations  (red  lines  in  Figure 
12a)  of  various  fits  of  the  measured  concentration  to  show  the  variability  generated  through  the  used  sim¬ 
ulation  approach. 


regularized  inversion  allows  us  to  evaluate  (1)  the  role  of 
(f>  in  the  transport  process,  (2)  whether  it  is  necessary  to 
include  (j)  variation  in  the  estimation  procedure,  (3)  the 
value  of  using  geophysical  data  to  develop  the  0  model, 
and  (4)  the  consequences  of  choices  made  in  parameteriz¬ 
ing  the  inversion  process  on  the  results  and  on  fitting  the 
measurements.  All  of  the  inversions  presented  here  were 
performed  using  seven  pilot  points  positioned  at  a  succession 
of  depths  at  a  given  x-y  coordinate  position  (shown  with 
black  dots  in  Figures  13-18).  A  layered  model  is  obtained  by 
kriging  with  very  large  lateral  correlation  lengths.  These 
choices  represent  a  trade-off  between  simplifying  the  prob¬ 
lem  and  increasing  the  number  of  pilot  points  and  problem 
dimensions,  where  choosing  the  more  complicated  approach 
would  (1)  significantly  increase  computing  demand,  (2) 
require  more  information  to  constrain/control  an  increasingly 
ill-posed  problem,  and  (3)  risk  estimation  of  a  model  with 
more  heterogeneity  than  is  needed  to  explain  the  data. 
Indeed,  the  large  hydrological  units  visible  at  the  BHRS 
from  other  sources  of  information  (Figure  3)  indicate  that  a 
simplified  approach  might  estimate  K  and  </>  structures  suffi¬ 
ciently  well  to  match  the  BTCs  to  a  first  order.  The  regulari¬ 
zation  for  the  inversion  process  was  imposed  on  the 
hydrological  model  using  second-derivative  smoothness 


constraints  in  the  x,  y,  and  z  directions.  Note  too,  that  we 
also  tested  a  Bayesian  approach  which  yielded  similar  results 
using  an  estimated  prior  mean  and  covariance  matrix. 

6.3.1.  Inversion  of  Hydraulic  Conductivity  Using  a 
Predefined  Constant  Effective  Porosity 

[42]  It  is  quite  common  to  set  to  a  constant  value  in  an 
inversion  procedure  to  estimate  K  because  the  influence  of 
the  <p  distribution  on  the  flow  and  transport  process  is  often 
assumed  to  be  small  compared  to  that  of  K.  Doing  this  also 
reduces  computing  demand.  Here  we  consider  two  cases 
where  (!)  is  assumed  known  and  (1)  equal  to  0.3,  and  then 
(2)  equal  to  0.2.  Both  values  can  be  considered  as  reasona¬ 
ble  estimates  for  the  aquifer  at  the  BHRS  based  on  the 
available  measurements  at  the  site. 

[43]  Figure  13  shows  the  resulting  K  and  the  predicted 
BTCs  when  (!)  is  set  to  0.3.  The  RMS  difference  between 
the  predicted  and  measured  BTCs  is  0.0305,  and  thus  can 
be  considered  as  relatively  good.  Also,  the  main  depth  var¬ 
iations  in  K  are  relatively  consistent  with  some  previous 
results  (Figures  9-11)  and  with  high-porosity  structures 
seen  in  </>  logs  (Figure  2)  below  838  m  between  wells  B3 
and  A1  and  imaged  by  geophysical  data  (Figure  3).  How¬ 
ever,  when  (j)  is  0.2  (Figure  14),  the  resulting  K  model 
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Figure  13.  Best  fit  model  obtained  through  a  regularized  least  square  inversion  to  estimate  K  distribu¬ 
tion:  (a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown 
in  Figure  5),  (b)  assumed  known  0,  and  (c)  estimated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e) 
simulated  concentration  after  9.3  and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 


yields  a  poorer  RMS  difference  (0.0484)  between  predicted 
and  measured  BTCs.  Also,  the  simulated  plume  crossing 
profile  B1-B4  at  various  times  (Figure  14)  is  not  qualita¬ 
tively  similar  to  the  time-lapse  GPR  tomography  (Figure  6). 
That  is,  a  difference  of  0. 1  in  the  assumed  known  0  has  a 
large  influence  on  the  results,  and  a  0  value  that  is  too  low, 
in  this  case  0.2,  does  not  allow  the  simulated  plume  to  sink 
enough  regardless  of  the  associated  K  distribution  estimated 


through  the  inversion  process.  Importantly,  these  results 
are  consistent  with  those  obtained  from  the  grid  search 
approach  applied  to  a  homogeneous  hydrological  model, 
where  low  </>  does  not  allow  good  fits  to  the  measured 
BTCs  (Figure  7).  Finally,  these  results  show  that  setting  0 
to  predefined  values,  which  is  quite  common  in  hydrologi¬ 
cal  studies,  can  significantly  influence  the  results  from  the 
inversion  process. 
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Figure  14.  Best  fit  model  obtained  through  a  regularized  least  square  inversion  to  estimate  K  distribu¬ 
tion:  (a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown 
in  Figure  5),  (b)  assumed  known  0,  and  (c)  estimated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e) 
simulated  concentration  after  9.3  and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 
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Figure  15.  Best  fit  model  obtained  through  a  regularized  least  square  inversion  to  estimate  K  and  6 
distributions  and  using  a  relatively  large  regularization  weighting  coefficient :  (a)  measured  (thin  lines) 
and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown  in  Figure  5),  (b)  estimated 
cj),  and  (c)  estimated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e)  simulated  concentration  after  9.3 
and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 


6.3.2.  Inversion  of  Hydraulic  Conductivity  and 
Effective  Porosity 

[44]  Deterministic  inversion  has  been  performed  to  esti¬ 
mate  both  K  and  0  for  the  following  two  cases.  First  (Figure 
15),  the  regularization  weighting  coefficient  (3  is  selected 
based  on  the  trade-off  between  data  fit  and  model  structure. 


This  means  the  structure,  or  more  exactly  the  values  at  the 
pilot  points,  are  enforced  for  continuity  (smoothness).  Sec¬ 
ond  (Figure  16),  a  smaller  regularization  coefficient  is 
applied  and  thus  the  variability  between  the  estimated  val¬ 
ues  at  the  pilot  points  may  be  higher.  Note  that  although 
kriging  controls  some  spatial  variability  between  the  pilot 
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Figure  16.  Best  fit  model  obtained  through  a  regularized  least  square  inversion  to  estimate  K  and  6 
distributions  and  using  a  relatively  low  regularization  weighting  coefficient:  (a)  measured  (thin  lines)  and 
calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown  in  Figure  5),  (b)  estimated  (b,  and 
(c)  estimated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e)  simulated  concentration  after  9.3  and  12 
d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 
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points,  the  variability  between  the  values  estimated  at  the 
pilot  points  is  controlled  only  by  the  regularization. 

[45]  For  the  first  case  (Figure  15),  the  RMS  difference 
(0.03 19)  between  the  estimated  and  measured  BTCs  is  rela¬ 
tively  good  and  the  K  and  cp  distributions  show  nice  agree¬ 
ment  with  previous  results  (Figures  9,  10,  11,  and  13)  and 
geophysical  data  (Figure  3)  and  <j)  logs  (Figure  2),  espe¬ 
cially  with  regard  to  increases  with  depth  at  ~838  m.  How- 
ever,  when  inverting  the  data  with  a  small  regularization 
weighting  parameter  (Figure  16),  the  estimated  K  and  tj> 
distributions  are  very  different.  In  contrast  to  Figure  15,  K 
and  <f>  distributions  in  Figure  16  show  strong  depth-related 
fluctuation  and  are  negatively  correlated.  Here  the  low  reg¬ 
ularization  parameter  has  enabled  more  variability  between 
the  values  estimated  at  the  pilot  points,  which  has  led  to  a 
completely  different  solution  with  a  lower  RMS  difference 
(0.0265).  However,  these  results  are  not  very  realistic  con¬ 
sidering  the  poor  similarity  in  variability  between  the  esti¬ 
mated  (j)  (Figure  16)  and  the  <j)  log  (Figure  2)  and  the 
geophysical  data  (Figures  3  and  6).  That  is,  the  choice  of  the 
regularization  parameter  may  result  in  very  different  solu¬ 
tions  even  when  structures  are  evaluated  at  a  relatively  large 
scale.  This  nonuniqueness  is  likely  driven  by  the  fact  that 
both  K  and  <p  distributions  are  estimated  through  the  inver¬ 
sion  process,  and  more  so  without  any  assumption  about  a 
relation  between  them.  So  this  again  shows  the  advantages 
of  and  need  to  constrain  the  estimation  of  hydrological  prop¬ 
erties  although,  at  the  same  time,  constraints  applied  on  the 
inversion  process  can  strongly  influence  the  result  and  thus 
need  to  be  reliable. 

6.3.3.  Inversion  of  Hydraulic  Conductivity  Using  a 
Predefined  3-D  Porosity  Distribution 

[46]  From  the  above  results  it  is  clear  that  estimation  of 
the  3-D  (f>  distribution  can  be  important  to  constrain  the 


inversion  of  hydrological  data  to  estimate  K.  Here  this  is 
done  by  using  a  <fi  realization  obtained  at  the  BHRS  from 
<p  logs  (e.g.,  Figure  2)  and  GPR  cross-hole  velocity  tomo¬ 
grams  with  a  simulated-annealing  approach  (e.g.,  Figure 
3b).  The  </>  realization  used  here  can  be  considered  as  repre¬ 
sentative  of  a  large  number  of  realizations  in  the  context  of 
this  study. 

[47]  Figure  1 7  shows  the  K  distribution  estimated  through 
the  regularized  inversion  when  a  simulation  of  <p  from  geo¬ 
physical  data  is  used  to  define  the  </>  model.  However,  the 
estimated  model  of  K  does  not  provide  a  satisfactory  RMS 
difference  (0.0759)  between  the  predicted  and  measured 
BTCs.  In  fact,  this  result  is  similar  to  that  obtained  when 
setting  the  <f>  equal  to  0.2  prior  to  the  inversion  (Figure  9) 
and  is  similar  to  the  results  from  the  grid  search  approach 
applied  to  a  homogeneous  model  (Figure  7).  In  all  these 
cases,  having  (j>  values  that  are  too  low  results  in  models 
that  poorly  match  the  measured  BTCs.  Given  these  results 
we  decided  to  reran  the  same  deterministic  inversion  after 
adding  0.1  to  the  estimated  <f>  distribution  based  on  the  pre¬ 
vious  grid  search  and  inversion  results.  In  this  case  (Figure 
1 8),  the  model  of  K  is  more  realistic  and  results  in  a  much 
smaller  RMS  difference  (0.0308)  between  predicted  and 
measured  BTCs.  Furthermore,  the  estimated  K  distribution 
now  agrees  well  with  previous  results  showing  similar  fea¬ 
tures  (Figures  9,  10,  13,  and  15)  and  shows  a  better  agree¬ 
ment  with  the  GPR  time-lapse  plume  images  (Figure  6) 
compared  to  Figure  17.  This  result  shows  that  the  structure 
estimated  from  geophysical  data  can  help  constrain  the  solu¬ 
tion  domain  and  is  not  meant  to  indicate  that  estimated  <p 
from  neutron  logs  or  other  geophysical  data  will  necessarily 
need  to  be  corrected  for  use  in  the  estimation  of  K.  In  fact, 
several  factors  may  explain  this  including  (1)  the  difficulty 
of  modeling  and  understanding  the  near-well  processes, 
especially  density-effects  during  the  injection  and  early 
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Figure  17.  Best  fit  model  obtained  through  a  regularized  least  square  inversion  to  estimate  K  distribution: 

(a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown  in  Figure  5), 

(b)  (p  from  geophysical  and  well  data,  and  (c)  estimated  K  distributions,  shown  along  profile  B6-B3 ;  (d-e) 
simulated  concentration  after  9.3  and  12  d  observed  along  profiles  A1-B3  and  B1-A1-B4,  respectively. 
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Figure  18.  Best  fit  model  obtained  through  a  regularized  least  square  inversion  to  estimate  K  distribu¬ 
tion:  (a)  measured  (thin  lines)  and  calculated  (thick  lines)  BTCs  in  20  depth  intervals  along  A1  (shown 
in  Figure  5),  (b)  0  from  geophysical  and  well  data  with  0.1  added,  and  (c)  estimated  K  distributions, 
shown  along  profile  B6-B3;  (d-e)  simulated  concentration  after  9.3  and  12  d  observed  along  profiles 
A1-B3  and  B1-A1-B4,  respectively. 


spreading  of  the  tracer,  (2)  the  importance  of  possible  high 
cf)  lenses  (some  visible  on  logs  between  838  and  834  m  in 
wells  B3  and  Al,  including  lenses  with  cj)  >  0.30)  that  may 
be  difficult  to  simulate  adequately  as  they  are  at  the  upper 
end  of  the  0  distribution  and  may  have  greater  occurrence 
than  shown  in  the  available  data  and  global  statistics  of  the 
site.  Uncertainties  in  these  areas  may  be  resolved  through 
future  density-driven  field,  laboratory,  and  numerical  stud¬ 
ies  and  through  coupling  the  hydrogeophysical  cj)  estimation 
and  the  tracer  data  inversion  together  to  ensure  that  a  signif¬ 
icant  near-well  structure  is  adequately  treated. 

7.  Summary  and  Conclusions 

[48]  The  objective  of  this  study  was  to  investigate  sev¬ 
eral  fundamental  issues  related  to  a  tracer  test  performed 
in  2001  at  the  BHRS  including  density  effects.  Our  interest 
was  to  improve  the  understanding  of  (1)  hydrologic  param¬ 
eters  influencing  such  transport  processes,  (2)  structural 
and  hydrological  parameter  information  that  can  be  gained 
from  tracer  concentration  and  head  data  and  from  other 
sources  of  geophysical  information  and  how  they  may  be 
consistent  or  complementary,  and  (3)  how  different  param¬ 
eter  estimation  methods  may  allow  various  results  from  the 
information  contained  in  the  data. 

[49]  With  respect  to  the  findings,  this  study  has  shown 
that,  at  the  BHRS  at  least,  0  has  an  important  role  in  the 
transport  and  sinking  of  solute.  In  fact,  through  the  estima¬ 
tion  methods  used  in  this  study  it  is  clear  that  choices  made 
about  how  the  (f>  distribution  is  defined  or  included  in  the 
estimation  approach  can  produce  strong  differences  in  the 
estimated  K  distribution.  In  addition,  if  0  is  defined  as  a  prior, 
some  values  preclude  finding  a  solution  with  acceptable  fit 


to  the  measurements.  This  finding  corroborates  results  from 
a  recent  synthetic  study  [Hu  et  al.,  2009],  showing  that 
variation  in  can  significantly  influence  solute  plume 
development. 

[50]  Furthermore,  although  hydrological  experiments, 
such  as  the  one  investigated  here,  are  known  to  be  relatively 
ill-posed  for  the  estimation  of  hydrological  property  distri¬ 
butions,  this  is  rarely  examined  in  the  literature.  This  study 
has  shown  how  choices  about  constraints  (assumptions, 
parameterization,  or  regularization),  which  are  commonly 
used  to  limit  effects  of  ill-posedness  in  the  estimation  proce¬ 
dure,  can  significantly  influence  the  estimated  hydrological 
parameters.  Also,  for  such  ill-posed  problems  where  non¬ 
uniqueness  can  be  expected,  the  comparison  of  multiple 
estimation  approaches  and  parameterizations  can  be  useful 
for  evaluating  which  inferred  hydrological  structures  may 
be  realistic,  and  the  extent  to  which  hydrological  structures 
with  different  scale  and/or  shape  may  allow  comparable  fits 
of  measured  BTCs. 

[51]  With  respect  to  the  use  of  geophysical  information  to 
improve  the  simulation  of  hydrological  processes,  here  simi¬ 
lar  hydrological  structures  seen  in  the  results  from  different 
parameter  estimation  approaches  are  in  good  agreement  with 
structures  in  cross-hole  GPR  velocity  tomograms  and  0  logs. 
Some  similarity  between  GPR  time-lapse  images  and  simu¬ 
lated  concentration  at  various  times  has  also  been  observed. 
Also,  this  study  has  shown  again  that  the  use  of  geophysical 
or  other  additional  information  is  a  key  to  potentially  reduc¬ 
ing  the  ill-posedness  present  in  the  hydrological  parameter 
estimation  process.  In  this  context,  the  structure  imaged  from 
geophysical  data  can  clearly  help  to  constrain  the  solution 
domain.  With  regard  to  constraining  the  tracer  data  inversion 
with  prior  estimated  <b  distribution,  the  strong  influence  of 
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the  hydrological  structure  and  processes  around  the  injection 
well  for  early  spreading  of  the  plume  require  further  field, 
laboratory,  and  numerical  investigations  to  improve  our 
understanding  of  such  processes  and  further  investigation  of 
how  to  ideally  couple  the  4>  estimation  process  with  the 
tracer  data  inversion. 

[52]  Finally,  the  analysis  of  this  tracer  test  data  may  bring 
some  new  insights  about  hydrological  unit  and  property  dis¬ 
tribution.  The  incomplete  fitting  of  BTCs,  which  in  fact,  is 
quite  common  to  some  degree  in  field  studies,  may  be 
because  of  difficulties  in  very  accurately  modeling  the  trans¬ 
port  process,  the  boundary  conditions,  and/or  the  details  of 
heterogeneity  of  hydrological  parameter  distributions.  Never¬ 
theless,  on  the  basis  of  the  inversion  results  here  and  on 
high-porosity  lenses  observed  in  $  logs,  the  presence  of  a 
relatively  higher  K  and  <j>  zone  below  838  m  compared  to 
above  838  m  seems  plausible.  At  this  point,  the  results 
from  multilevel  slug  tests  performed  very  recently  [ Cardiff 
et  al.,  2011],  show  consistently  higher  K  between  836  and 
838  m  than  between  838  m  and  839.5  m,  although  the 
obtained  K  values  are  higher  than  observed  in  this  study. 
The  same  study,  Cardiff  el  al.  [2011]  shows  that  the  com¬ 
monly  used  assumption  of  a  strong  correlation  between  <j> 
and  K  in  unconsolidated  coarse  sedimentary  aquifers  may 
not  be  fully  applicable  at  the  BHRS.  This,  as  well  as  com¬ 
parison  and  interpretation  of  the  various  hydrogeological 
structures  recognized  at  the  BHRS  will  be  the  subject  of 
future  research.  In  a  more  general  context,  insights  into  the 
use  of  the  methods  examined  here  can  hopefully  help  future 
tracer  test  experiments  with  choices  in  acquisition  and 
processing,  strategies  for  multiple  geophysical  and  hydro- 
logical  measurement  acquisitions,  and  choices  made  during 
modeling  and  inversion  of  the  experimental  data. 
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